home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / answrbok / 6_10.lha / 6_10 / 6_10div.c < prev    next >
Text File  |  1993-08-08  |  14KB  |  534 lines
  1. * Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
  2. * The C++ Answer Book */
  3. * Tony Hansen */
  4. * All rights reserved. */
  5. *
  6.    Divide u[1..m+n] by v[1..n] to
  7.    form q[0..m] and r[1..n]
  8.  
  9.    Based on:
  10.  
  11.    The Art of Computer Programming, volume 2
  12.    D. Knuth, Section 4.3.1, Algorithm D and
  13.    exercise 16.
  14. /
  15. include <lint.h>
  16. include <debug.h>    /* DELETE */
  17.  
  18. oid dodivmod(LINT u, LINT v,
  19.    LINT "ient, LINT &remainder)
  20.  
  21.    int negquotient = 0, negremainder = 0;
  22.  
  23.                                                             { if (debug > 1) cerr << "\n"; /* DELETE */ }
  24.                                                             { if (debug > 1) cerr << form("u=%4.4x%4.4x%4.4x%4.4x", u.s[0], u.s[1], u.s[2], u.s[3]) << "\n"; /* DELETE */ }
  25.                                                             { if (debug > 1) cerr << form("v=%4.4x%4.4x%4.4x%4.4x", v.s[0], v.s[1], v.s[2], v.s[3]) << "\n"; /* DELETE */ }
  26.    /*
  27. Make u (the dividend) positive.
  28. The sign of the remainder will
  29. match the sign of the dividend.
  30.    */
  31.    if (u.isneg())
  32. {
  33. negremainder = 1;
  34. u = -u;
  35.  if (debug) cerr << "PT 1\n";    /* DELETE */ }
  36.                                                             { if (debug > 1) cerr << "neg remainder\n"; /* DELETE */ }
  37. }
  38.    else    /* DELETE */
  39. {    /* DELETE */
  40.  if (debug) cerr << "PT 2\n";    /* DELETE */ }
  41. }    /* DELETE */
  42.  
  43.    /*
  44. Make v (the divisor) positive.
  45. The sign of the quotient will be
  46. negative if the sign of the divisor
  47. and dividend do not match, else
  48. positive.
  49.    */
  50.    if (v.isneg())
  51. {
  52.  if (debug) cerr << "PT 3\n";    /* DELETE */ }
  53. negquotient = !negremainder;
  54. v = -v;
  55. }
  56.    else
  57. {    /* DELETE */
  58.  if (debug) cerr << "PT 4\n";    /* DELETE */ }
  59. negquotient = negremainder;
  60. }    /* DELETE */
  61.                                                             { if (debug > 1) cerr << (negquotient ? "negquotient\n" : ""); /* DELETE */ }
  62.  
  63.    // set local variables
  64.    for (int m_n = 4, uoffset = 0;
  65.  uoffset < 4 && u.s[uoffset] == 0;
  66.  m_n--, uoffset++)
  67. {    /* DELETE */
  68.  if (debug) cerr << "PT 5\n";    /* DELETE */ }
  69. ;
  70. }    /* DELETE */
  71.  
  72.    for (int n = 4, voffset = 0;
  73.  voffset < 4 && v.s[voffset] == 0;
  74.  n--, voffset++)
  75. {    /* DELETE */
  76.  if (debug) cerr << "PT 6\n";    /* DELETE */ }
  77. ;
  78. }    /* DELETE */
  79.  
  80.    int m = m_n - n;
  81.                                                             { if (debug > 1) cerr << "m=" << m <<", n=" << n << ", m_n=" << m_n << "\n"; /* DELETE */ }
  82.  
  83.    /*
  84. If n == 0, then division by zero
  85.    */
  86.    if (n == 0)
  87. {
  88.  if (debug) cerr << "PT 7\n";    /* DELETE */ }
  89.                                                             { if (debug > 1) cerr << "division by zero\n"; /* DELETE */ }
  90. error("division by zero!");
  91. quotient = remainder = u;
  92. return;
  93. }
  94.  
  95.    /*
  96. For n == 1, use simpler algorithm
  97.    */
  98.    else if (n == 1)
  99. {
  100.  if (debug) cerr << "PT 8\n";    /* DELETE */ }
  101.                                                             { if (debug > 1) cerr << "n==1\n"; /* DELETE */ }
  102. // See Exercise 16 after Section 4.3.1
  103. LINT q = 0;
  104. LINT_Ltype prevu = 0;
  105. LINT_type v1 = v.s[3];
  106.  
  107. for (int r = uoffset; r < m_n + uoffset; r++)
  108.     {
  109.  if (debug) cerr << "PT 9\n";    /* DELETE */ }
  110.     LINT_Ltype t = u.s[r] + prevu * LINT_base;
  111.     LINT_type tmpq = LINT_type (t / v1);
  112.     q.s[r] = tmpq;    // % LINT_base
  113.     prevu = t - v1 * tmpq;
  114.     }
  115.  
  116. if (negquotient)
  117.     {    /* DELETE */
  118.  if (debug) cerr << "PT 10\n";    /* DELETE */ }
  119.     quotient = -q;
  120.     }    /* DELETE */
  121. else
  122.     {    /* DELETE */
  123.  if (debug) cerr << "PT 11\n";    /* DELETE */ }
  124.     quotient = q;
  125.     }    /* DELETE */
  126.  
  127. if (negremainder)
  128.     {    /* DELETE */
  129.  if (debug) cerr << "PT 12\n";    /* DELETE */ }
  130.     remainder = -prevu;
  131.     }    /* DELETE */
  132. else
  133.     {    /* DELETE */
  134.  if (debug) cerr << "PT 13\n";    /* DELETE */ }
  135.     remainder = prevu;
  136.     }    /* DELETE */
  137. return;
  138. }
  139.  
  140.    /*
  141. Degenerate case of length(u) < length(v)
  142. i.e., m < 0, implying that u < v
  143.    */
  144.    else if (m_n < n)
  145. {
  146.  if (debug) cerr << "PT 14\n";    /* DELETE */ }
  147.                                                             { if (debug > 1) cerr << "m_n < n\n"; /* DELETE */ }
  148. quotient = 0L;
  149. if (negremainder)
  150.     {    /* DELETE */
  151.  if (debug) cerr << "PT 15\n";    /* DELETE */ }
  152.     remainder = -u;
  153.     }    /* DELETE */
  154. else
  155.     {    /* DELETE */
  156.  if (debug) cerr << "PT 16\n";    /* DELETE */ }
  157.     remainder = u;
  158.     }    /* DELETE */
  159. return;
  160. }
  161.  
  162.    /*
  163. Degenerate case of length(u) == length(v)
  164. i.e., m == 0, possibly implying that
  165. u < v or u == v
  166.    */
  167.    else if (m_n == n)
  168. {
  169.  if (debug) cerr << "PT 17\n";    /* DELETE */ }
  170.                                                             { if (debug > 1) cerr << "m_n == n\n"; /* DELETE */ }
  171. int cmp = LINT_cmp(&u.s[uoffset], &v.s[voffset], m_n);
  172. if (cmp < 0)    // u < v
  173.     {
  174.  if (debug) cerr << "PT 18\n";    /* DELETE */ }
  175.                                                             { if (debug > 1) cerr << "cmp < 0\n"; /* DELETE */ }
  176.     quotient = 0L;
  177.     if (negremainder)
  178.     {    /* DELETE */
  179.  if (debug) cerr << "PT 19\n";    /* DELETE */ }
  180.     remainder = -u;
  181.     }    /* DELETE */
  182.     else
  183.     {    /* DELETE */
  184.  if (debug) cerr << "PT 20\n";    /* DELETE */ }
  185.     remainder = u;
  186.     }    /* DELETE */
  187.     return;
  188.     }
  189.  
  190. else if (cmp == 0)    // u == v
  191.     {
  192.  if (debug) cerr << "PT 21\n";    /* DELETE */ }
  193.                                                             { if (debug > 1) cerr << "cmp == 0\n"; /* DELETE */ }
  194.     if (negquotient)
  195.     {    /* DELETE */
  196.  if (debug) cerr << "PT 22\n";    /* DELETE */ }
  197.     quotient = -1L;
  198.     }    /* DELETE */
  199.     else
  200.     {    /* DELETE */
  201.  if (debug) cerr << "PT 23\n";    /* DELETE */ }
  202.     quotient = 1L;
  203.     }    /* DELETE */
  204.     remainder = 0L;
  205.     return;
  206.     }
  207. else    /* DELETE */
  208.     {    /* DELETE */
  209.  if (debug) cerr << "PT 24\n";    /* DELETE */ }
  210.     }    /* DELETE */
  211. }
  212.    else    /* DELETE */
  213. {    /* DELETE */
  214.  if (debug) cerr << "PT 25\n";    /* DELETE */ }
  215. }    /* DELETE */
  216.  
  217.    /*
  218. Now call out all of the guns from Knuth
  219.    */
  220.  
  221.    /*
  222. D1(a) [Normalize.]
  223.     Set d <- b/(v1+1).
  224.    */
  225.    LINT_type d = LINT_base / (v.s[voffset] + 1);
  226.                                                             { if (debug > 1) cerr << "d=" << form("0x%4.4x", d) << "\n"; /* DELETE */ }
  227.  
  228.    /*
  229. D1(b) [Normalize.]
  230.     Set (u[0]u[1]...u[m+n]) to
  231.     (u[1]u[2]...u[m+n])*d
  232.    */
  233.    LINT_type uv[5];
  234.    LINT_type k;
  235.    if (d == 1)
  236. {
  237.  if (debug) cerr << "PT 26\n";    /* DELETE */ }
  238.                                                             { if (debug > 1) cerr << "d == 1, copy old u value\n"; /* DELETE */ }
  239. // copy old value
  240. uv[0] = 0;
  241. memcpy((char*)&uv[1], (char*)&u.s[uoffset],
  242.     m_n * sizeof(LINT_type));
  243. }
  244.  
  245.    else
  246. {
  247.  if (debug) cerr << "PT 27\n";    /* DELETE */ }
  248.                                                             { if (debug > 1) cerr << "multiply u by d\n"; /* DELETE */ }
  249. // multiply u by d
  250. k = 0;
  251. for (int Oi = m_n - 1 + uoffset, i = m_n;
  252.      i > 0; Oi--, i--)
  253.     {
  254.  if (debug) cerr << "PT 28\n";    /* DELETE */ }
  255.     LINT_Ltype t = u.s[Oi] * d + k;
  256.     uv[i] = t;        // % LINT_base
  257.     k = t / LINT_base;
  258.                                                             { if (debug > 1) cerr << "done setting u.s[" << Oi << "], uv[" << i << "] = " << form("0x%4.4x", uv[i]) << "\n"; /* DELETE */ }
  259.     }
  260. uv[i] = k;
  261. }
  262.                                                             { if (debug > 1) { cerr << "uv=0x"; for (int jj=0; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
  263.  
  264.    /*
  265. D1(c) [Normalize.]
  266.     Set (v[1]v[2]...v[n]) to
  267.     (v[1]v[2]...v[n]) * d
  268.    */
  269.    LINT_type vv[4];
  270.    if (d == 1)
  271. {
  272.  if (debug) cerr << "PT 29\n";    /* DELETE */ }
  273.                                                             { if (debug > 1) cerr << "d == 1, copy old v value\n"; /* DELETE */ }
  274. // copy old value
  275. memcpy((char*)&vv[0], (char*)&v.s[voffset],
  276.     n * sizeof(LINT_type));
  277. }
  278.  
  279.    else
  280. {
  281.  if (debug) cerr << "PT 30\n";    /* DELETE */ }
  282.                                                             { if (debug > 1) cerr << "multiply v by d\n"; /* DELETE */ }
  283. // multiply v by d
  284. k = 0;
  285. for (int Oi = n - 1 + voffset, i = n - 1;
  286.      i >= 0; Oi--, i--)
  287.     {
  288.  if (debug) cerr << "PT 31\n";    /* DELETE */ }
  289.     LINT_Ltype t = v.s[Oi] * d + k;
  290.     vv[i] = t;        // % LINT_base;
  291.     k = t / LINT_base;
  292.                                                             { if (debug > 1) cerr << "done setting v.s[" << Oi << "], vv[" << i << "] = " << form("0x%4.4x", vv[i]) << "\n"; /* DELETE */ }
  293.     }
  294. }
  295.                                                             { if (debug > 1) { cerr << "vv=0x"; for (int jj=0; jj < n; jj++) cerr << form("%4.4x", vv[jj]); cerr << "\n"; } /* DELETE */ }
  296.  
  297.    /*
  298. D2 [Initialize j]
  299.     Set j <- 0
  300.  
  301. D7 [Loop on j]
  302.     Set j <- j+1
  303.     Loop if j <= m
  304.    */
  305.    LINT_type v1 = vv[0];
  306.    LINT_type v2 = vv[1];
  307.                                                             { if (debug > 1) cerr << "v1=" << form("0x%4.4x", v1) << ", v2=" << form("0x%4.4x", v2) << "\n"; /* DELETE */ }
  308.    LINT_type qv[5];
  309.  
  310.    for (int j = 0; j <= m; j++)
  311. {
  312.  if (debug) cerr << "PT 32\n";    /* DELETE */ }
  313.                                                             { if (debug > 1) cerr << "j=" << j << "\n"; /* DELETE */ }
  314. /*
  315.     D3 [Calculate q^]
  316.     If u[j] == v[1]
  317.         Set q^ <- base - 1
  318.     Else
  319.         Set q^ <-
  320.         (u[j] * base + u[j+1]) / v[1]
  321.     While v[2] * q^ >
  322.        (u[j] * base + u[j+1] - q^ * v[1]) *
  323.         base + u[j+2]
  324.         Set q^ <- q^ - 1
  325. */
  326. LINT_Ltype q_hat;
  327.  
  328. if (uv[j] == v1)
  329.     {    /* DELETE */
  330.  if (debug) cerr << "PT 33\n";    /* DELETE */ }
  331.     q_hat = LINT_base - 1;
  332.                                                             { if (debug > 1) cerr << "uv[" << j << "]=" << form("0x%4.4x", uv[j]) << " == v1, q^=" << form("0x%4.4x", q_hat) << "\n"; /* DELETE */ }
  333.     }    /* DELETE */
  334. else
  335.     {    /* DELETE */
  336.  if (debug) cerr << "PT 34\n";    /* DELETE */ }
  337.     q_hat = (uv[j] * LINT_base + uv[j+1]) / v1;
  338.                                                             { if (debug > 1) cerr << "uv[" << j << "]=" << form("0x%4.4x", uv[j]) << " != v1, q^=" << form("0x%4.4x", q_hat) << "\n"; /* DELETE */ }
  339.     }    /* DELETE */
  340.  
  341. LINT_Ltype u_j = uv[j];
  342. LINT_Ltype u_j1 = uv[j+1];
  343. LINT_Ltype u_j2 = uv[j+2];
  344.  
  345. for ( ; ; q_hat--)
  346.     {
  347.  if (debug > 1) cerr << "q^=" << form("0x%4.4x", q_hat) << "\n"; /* DELETE */ }
  348.     /*
  349.     The following statements build up the
  350.     following test which combines with the
  351.     for loop above to give the while
  352.     statement mentioned in the last comment.
  353.     if ((v2 * q_hat) <=
  354.         ((u_j * LINT_base + u_j1 - v1 *
  355.           q_hat) * LINT_base + u_j2))
  356.         q^--;
  357.     */
  358.  
  359.     // u_j_q_hat = u_j * LINT_base + u_j1
  360.     union {
  361.     LINT_Ltype L[2];
  362.     LINT_type S[4];
  363.     } u_j_q_hat;
  364.     u_j_q_hat.L[0] = 0;
  365.     u_j_q_hat.S[2] = u_j;
  366.     u_j_q_hat.S[3] = u_j1;
  367.  
  368.     // ... - v1 * q_hat
  369.     u_j_q_hat.L[1] -= v1 * q_hat;
  370.     if (u_j_q_hat.S[2] != 0)
  371.     break;
  372.  
  373.     // ... * LINT_base
  374.     u_j_q_hat.S[2] = u_j_q_hat.S[3];
  375.  
  376.     // ... + u_j2
  377.     u_j_q_hat.S[2] = u_j2;
  378.     if ((v2 * q_hat) <= u_j_q_hat.L[1])
  379.     break;
  380.  if (debug) cerr << "PT 35\n";    /* DELETE */ }
  381.  if (debug > 1) cerr << "q^--\n"; /* DELETE */ }
  382.     }
  383.  
  384. /*
  385.     D4 [Multiply and subtract.]
  386.     Replace u[j..j+n] by
  387.         u[j..j+n] -
  388.          q^ * v[1..n]
  389. */
  390. // set nv <- q^ * (v[1..n])
  391. LINT_type nv[5];
  392. k = 0;
  393. for (int dl = n, vl = n-1; vl >= 0; vl--, dl--)
  394.     {
  395.  if (debug) cerr << "PT 36\n";    /* DELETE */ }
  396.     LINT_Ltype t = vv[vl] * q_hat + k;
  397.     nv[dl] = t;        // % LINT_base;
  398.     k = t / LINT_base;
  399.                                                             { if (debug > 1) cerr << "done setting nv[" << dl << "] = " << form("0x%4.4x", nv[dl]) << ", vv[" << vl << "] = " << form("0x%4.4x", vv[vl]) << "\n"; /* DELETE */ }
  400.     }
  401. nv[0] = k;
  402.                                                             { if (debug > 1) cerr << "done setting nv[" << dl << "] = " << form("0x%4.4x", nv[dl]) << "\n"; /* DELETE */ }
  403.                                                             { if (debug > 1) { cerr << "nv=0x"; for (int jj=0; jj <= n; jj++) cerr << form("%4.4x", nv[jj]); cerr << "\n"; } /* DELETE */ }
  404.  
  405. // subtract nv[0..n] from u[j..j+n]
  406. int borrow = 0;
  407. int ul = j + n;
  408.                                                             { if (debug > 1) cerr << "ul=" << ul << ", dl=" << n << "\n"; /* DELETE */ }
  409. for (dl = n; dl >= 0; dl--, ul--)
  410.     {
  411.  if (debug) cerr << "PT 37\n";    /* DELETE */ }
  412.     LINT_Ltype t = uv[ul] - nv[dl] - borrow;
  413.     uv[ul] = t;        // % LINT_base
  414.     borrow = (t / LINT_base) ? 1 : 0;
  415.                                                             { if (debug > 1) cerr << "done setting uv[" << ul << "] = " << form("0x%4.4x", uv[ul]) << "\n"; /* DELETE */ }
  416.     }
  417.                                                             { if (debug > 1) { cerr << "uv=0x"; for (int jj=0; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
  418.                                                             { if (debug > 1) cerr << "borrow=" << borrow << "\n"; /* DELETE */ }
  419.  
  420. /*
  421.     D5 [Test remainder]
  422.     Set q[j] <- q^
  423.     if u[j] < 0
  424.  
  425.     D6 [Add back]
  426.         add 0v[1..n] back to u[j..j+n]
  427.         q[j]--
  428. */
  429. qv[j] = q_hat;
  430.                                                             { if (debug > 1) cerr << "qv[" << j << "] = " << form("0x%4.4x", qv[j]) << "\n"; /* DELETE */ }
  431. if (borrow != 0)
  432.     {
  433.  if (debug) cerr << "PT 38\n";    /* DELETE */ }
  434.                                                             { if (debug > 1) cerr << "borrow != 0, add v back to u\n"; /* DELETE */ }
  435.     for (k = 0, ul = j + n, vl = n - 1;
  436.      vl >= 0;
  437.      vl--, ul--)
  438.     {
  439.  if (debug) cerr << "PT 39\n";    /* DELETE */ }
  440.     LINT_Ltype t = uv[ul] + vv[vl] + k;
  441.     uv[ul] = t;    // % LINT_base
  442.                                                             { if (debug > 1) cerr << "done setting uv[" << ul << "] = " << form("0x%4.4x", uv[ul]) << "\n"; /* DELETE */ }
  443.     k = t / LINT_base;
  444.     }
  445.  
  446.     uv[j] += k;
  447.                                                             { if (debug > 1) cerr << "done setting uv[" << ul << "] = " << form("0x%4.4x", uv[ul]) << "\n"; /* DELETE */ }
  448.                                                             { if (debug > 1) { cerr << "uv=0x"; for (int jj=0; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
  449.     qv[j]--;
  450.                                                             { if (debug > 1) cerr << "qv[" << j << "] = " << form("0x%4.4x", qv[j]) << "\n"; /* DELETE */ }
  451.     }
  452. else    /* DELETE */
  453.     {    /* DELETE */
  454.  if (debug) cerr << "PT 40\n";    /* DELETE */ }
  455.     }    /* DELETE */
  456. }
  457.  
  458.    /*
  459. D8 [Unnormalize]
  460.     q[0..m] is quotient
  461.     u[m+1..m+n] / d is remainder
  462.    */
  463.    LINT q = 0;
  464.    memcpy((char*)&q.s[3-m], (char*)&qv[0],
  465. (m+1) * sizeof(LINT_type));
  466.                                                             { if (debug > 1) cerr << "q = " << q << "\n"; /* DELETE */ }
  467.    if (negquotient)
  468. {    /* DELETE */
  469.  if (debug) cerr << "PT 41\n";    /* DELETE */ }
  470. quotient = -q;
  471. }    /* DELETE */
  472.  
  473.    else
  474. {    /* DELETE */
  475.  if (debug) cerr << "PT 42\n";    /* DELETE */ }
  476. quotient = q;
  477. }    /* DELETE */
  478.  
  479.    // divide u[m+1..m+n] by d
  480.    LINT r = 0;
  481.                                                             { if (debug > 1) { cerr << "uv[m+1..m+n]=0x"; for (int jj=m+1; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
  482.    if (d == 1)        // nothing special to do
  483. {    /* DELETE */
  484.  if (debug) cerr << "PT 43\n";    /* DELETE */ }
  485. memcpy((char*)&r.s[4 - n], (char*)&uv[m + 1],
  486.     n * sizeof(LINT_type));
  487. }    /* DELETE */
  488.  
  489.    else
  490. {
  491.  if (debug) cerr << "PT 44\n";    /* DELETE */ }
  492. LINT_Ltype prevu = 0;
  493.  
  494. // do division by single digit
  495. for (int rl = 4 - n, ul = m + 1;
  496.      ul <= m_n;
  497.      ul++, rl++)
  498.     {
  499.  if (debug) cerr << "PT 45\n";    /* DELETE */ }
  500.     LINT_Ltype t = uv[ul] + prevu * LINT_base;
  501.     LINT_type tmpq = t / d;
  502.     r.s[rl] = tmpq;    // % LINT_base
  503.     prevu = t - d * tmpq;
  504.     }
  505. }
  506.  
  507.                                                             { if (debug > 1) cerr << "r = " << r << "\n"; /* DELETE */ }
  508.    if (negremainder)
  509. {    /* DELETE */
  510.  if (debug) cerr << "PT 46\n";    /* DELETE */ }
  511. remainder = -r;
  512. }    /* DELETE */
  513.  
  514.    else
  515. {    /* DELETE */
  516.  if (debug) cerr << "PT 47\n";    /* DELETE */ }
  517. remainder = r;
  518. }    /* DELETE */
  519.  
  520.  
  521. INT operator%(LINT u, LINT v)
  522.  
  523.    LINT quot, rem;
  524.    dodivmod(u, v, quot, rem);
  525.    return rem;
  526.  
  527.  
  528. INT operator/(LINT u, LINT v)
  529.  
  530.    LINT quot, rem;
  531.    dodivmod(u, v, quot, rem);
  532.    return quot;
  533.  
  534.